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Radiant  heat  transfer  plays  an  important  role  in  the  distribution  of  cell  temperature  and  current  density 
in  solid  oxide  fuel  cells  (SOFC).  The  objective  of  this  paper  is  to  introduce  a  mathematical  model  of  view 
factors  for  radiation  heat  exchange  in  an  in-house  longitudinally  distributed  SOFC  model.  A  differential 
view  factor  model  is  first  developed  for  planar  and  tubular  SOFC  configurations,  but  is  found  invalid  when 
the  infinitesimal  element  size  is  comparable  to  the  characteristic  size.  Then,  a  finite-difference  view  factor 
model  is  developed  to  solve  the  problem  of  discontinuities  in  the  differential  view  factor  model.  Starting 
from  a  classical  problem  of  convective  and  radiant  heat  transfer  for  a  transparent  gas  flow  in  a  gray-wall 
tube,  a  fast  and  accurate  computation  is  available  for  the  finite-difference  view  factor  model  without 
extra  mathematical  derivations  of  the  governing  equations.  Compared  to  the  simple  modeling  which  only 
takes  into  account  the  surface-to-surface  radiation  exchange  between  two  directly  opposed  elements, 
the  detailed  radiation  model  based  on  analytical  view  factors  predicts  more  uniform  distribution  of  cell 
temperature  and  current  density  in  the  overall  SOFC  modeling. 

©  2010  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

Solid  oxide  fuel  cells  (SOFCs)  are  a  type  of  high-temperature 
fuel  cells.  Planar  and  tubular  configurations  are  the  two  most  com¬ 
mon  SOFC  designs.  With  the  advantages  of  short  current  path  and 
lower  component  fabrication  cost,  the  planar  design  is  usually  the 
preferred  design,  but  it  suffers  from  sealing  problems  caused  by 
thermal  expansion  at  high-temperature.  Until  now,  planar  design 
has  been  widely  used  in  anode-supported  intermediate  temper¬ 
ature  SOFCs  [1].  The  tubular  design  is  superior  with  respect  to 
sealing,  but  often  has  to  operate  in  the  high  temperature  range 
of  900-1 000  °C  to  maintain  acceptable  ohmic  polarization  due  to 
long  current  paths.  Some  other  advanced  configurations,  such  as 
the  Siemens  Westinghouse  design  of  flattened  tubular  SOFC  and  the 
integrated  planar  SOFC  pioneered  by  Rolls-Royce  (which  combines 
the  advantages  of  the  planar  and  tubular  designs  by  introducing 
the  multi-cell  membrane  electrode  assembly  (MEA)  concept  [2]) 
can  be  considered  as  a  variation  of  the  two  basic  planar  and  tubular 
designs  [1]. 
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Radiation  heat  transfer  plays  an  important  role  in  the  distribu¬ 
tion  of  cell  temperature  and  current  density  in  SOFCs.  A  SOFC  model 
with  detailed  radiation  model  is  necessary  for  accurate  predic¬ 
tion  of  cell  performance.  A  spectral  radiative  heat  transfer  analysis 
within  a  computational  fluid  dynamics  (CFD)  SOFC  model  is  pow¬ 
erful  but  usually  computationally  expensive  [3].  A  simple  method 
of  modeling  the  surface-to-surface  radiation  exchange  is  generally 
acceptable  in  thermal  models  of  SOFCs.  However,  most  of  state-of- 
the-art  radiation  models  just  considered  the  radiation  exchange 
between  two  directly  opposed  differential  elements,  where  the 
view  factor  between  the  two  finite  elements  is  considered  identi¬ 
cal  to  that  between  two  infinite  diffuse  interchange  surfaces  [4-6]. 
Some  other  papers  includes  radiation  exchange  among  three  adja¬ 
cent  control  volumes,  but  the  constant  view  factors  are  strongly 
dependent  on  the  number  of  discrete  grids  [7,8].  Aiming  at  system- 
level  analysis,  an  in-house  multi-level  simulation  platform  for  solid 
oxide  fuel  cell  and  gas  turbine  hybrid  generation  system  was 
recently  developed  in  gPROMS,  a  commercial  advanced  process 
modeling  software  [9].  The  objective  of  this  paper  is  to  introduce 
an  analytical  model  of  view  factors  in  our  longitudinally  distributed 
models  of  planar  and  tubular  SOFCs  [10]. 

Fig.  1  shows  the  schematic  diagram  of  a  section  of  planar  and 
tubular  SOFCs.  For  planar  SOFCs  (PSOFC),  Wch  and  Dch  are  the  width 
and  depth  of  the  rectangular  flow  channel  enclosed  between  the 
positive  electrolyte  negative  (PEN)  and  the  interconnect  (CON).  The 
geometry  of  tubular  SOFCs  (TSOFC)  is  that  of  the  commercial  West- 
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Fig.  1.  Schematic  diagram  of  a  section  of  planar  and  tubular  SOFCs. 
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Fig.  2.  Configuration  for  diffuse  interchange  in  planar  SOFC. 


inghouse  cell  [1].  Fuel  flows  along  outside  of  the  tube,  while  air  is 
first  preheated  via  the  air  supply  tube  (AST)  and  then  flows  back 
into  the  annular  area  between  the  cathode  and  the  air  supply  tube. 
In  TSOFC,  rAST  in  and  rAST  out  are  the  inner  and  outer  radius  of  the 
AST,  rin  and  rout  are  the  inner  and  outer  radius  of  the  cell  tube.  The 
analytical  view  factors  between  the  solid  phases,  i.e.  PEN  and  CON 
in  PSOFCs,  and  PEN  and  AST  in  TSOFCs  will  be  discussed  in  the  next 
section. 


Define  X=Dch/Wch,  Z=z/Wch,  Y=  |Z2-Z!  |, 

92Fi  3  1 

dFdPENl ,  dCON2  =  dfdPENl ,  ds2  =  PCh  dz2  =  -  —f\  (X,  Y)dZ2  ( 1 ) 

OZ\  OZ2  7T 


dFdCONl ,  dPEN2  =  PEN’2  dFdPEN2;  dCONl  =  - 


dACoN,i 


tt(  1  +  2X)' 


7i  (X,  Y)dZ2 


(2) 


2.  Analytical  differential  view  factors 

Because  the  cell  length  (L)  is  generally  much  larger  than  the 
size  in  other  directions,  only  the  temperature  distribution  along 
the  axial  coordinate  ze[0,L]  is  considered  in  our  longitudinally  dis¬ 
tributed  models  of  both  planar  and  tubular  SOFCs. 


2.1.  PSOFC  configuration 


dFdCON  1 ,  dCON2 


WchPch 

Wch  +  2  Dch 


^1,2  2  a2Fi,3  \ 
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dz2 


where 
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2  n  (1  +Y2)(X2  +  Y2) 


The  configuration  for  diffusive  interchange  in  PSOFC  is  shown  in 
Fig.  2.  The  view  factor  between  two  identical,  parallel,  and  directly 
opposed  finite  rectangles  (Fi>2),  and  the  view  factor  between  two 
finite  rectangles  of  the  same  length,  having  one  common  edge,  and 
90°  from  each  other  (F13),  can  be  easily  found  in  the  view  factor 
catalogue  [11].  For  the  non-concave  PEN  surface,  dFdPEN1  dPEN2  =  0. 
According  to  the  reciprocity  and  additivity  rule,  the  view  factor 
between  two  infinitesimal  elements  can  be  directly  calculated  from 
the  second-order  differentiation  of  the  finite  view  factors,  Fii2  and 
Fii3,  as  shown  below. 
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Fig.  3.  Configuration  for  diffuse  interchange  in  tubular  SOFC. 


In  all  expressions  in  which  inverse  trigonometric  functions 
appear,  the  principal  value  is  to  be  taken,  i.e.  for  any  argument  x, 
— 7r/2  <  tg_1x  <  7t/2.  The  view  factor  from  the  infinitesimal  element 
to  the  two  opening  ends  of  the  groove  (Yi  =z/Wch,  Y2  =  (L-z)/Wch) 
is  related  to  the  first-order  differentiation  of  F\  2  and  Fii3: 


^dPEN,  end 
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y2  =  -rsin0m,  dy2  =  -rirdr/r2/(r2  -  r2)1/2,  re[r,  r+dr].  Line  41 
is  symmetric  to  line  23  along  the  x  axis,  and  arc  ^  34  is  in  the 
opposite  direction  of  arc  ^  21,  as  the  radius  is  equal  to  r  +  dr. 

When  defining  D  =  rin/rAST)0Ut,  Z=z/rAST,  out,  |Z2-Zi  |, 

A  =  (D  + 1  )2  +  H2,  B  =  (D  -  1  )2  +  Ft2,  the  view  factor  FdPENend  between 
the  ring  element  on  the  interior  of  the  PEN  cylinder  and  the  annular 
end  of  channel  (z2  =  0,  Ff=Z)  is 

rr=r2 

^dPEN,  end  =  /  dFdl*  d2 

J  r=r-[ 


1 

7 xD 


2 D2+H2  ,  a/(4D2  +H2)(D2  -  1) 
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D+T  tan  V  B(D  +  1) 


tan-’ 

H2  +  2(D2  —  1) 


(11) 


Similarly,  the  view  factor  FdASToutend  between  the  ring  element 
on  the  exterior  of  the  AST  cylinder  and  the  annular  end  of  the 
channel  (z2  =  0,  H=Z)  is 


^ dAST,  out,  end  — 


H  •  tan 


D-  1  ,  Jd2-\ 

+  tan-1  v 


D  +  l 


H 
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2.2.  TSOFC  configuration 


As  shown  in  Fig.  3,  for  the  configuration  of  coaxial  cylinders, 
the  view  factor  dFdl*  d2  from  the  differential  element  (dAi*)  at  the 
top  end  of  the  interior  surface  of  the  outer  cylinder  to  the  differen¬ 
tial  annular  element  (dA2)  at  the  base  of  the  outer  cylinder  can  be 
calculated  by  the  contour-integral  method  [12].  In  the  cylindrical 
coordinate  system, 


dFdU'd2=  2^ 


(Y2  -y\)dz2  ~(Z2  -*i  W2 

-12+23+^34+41  (*2  -  x\  f  +  (y2  -  Yl  )2  +  (Z2  ~  ^1  )2 

(10) 


The  view  factor  between  dA!  *  and  the  exterior  surface  of  the 
inner  cylinder  (A3)  and  the  interior  surface  of  the  outer  cylin¬ 
der  {A-[ ),  Fdl*  3  and  Fdri  can  also  be  obtained  by  calculating  the 
line  integral  along  aa'b'c'cba  and  deffe'd'd,  as  shown  in  Fig.  3.  The 
expression  of  Fdl*3  can  be  found  in  the  view  factor  catalogue  [11  ]. 
Then,  the  view  factors  between  the  ring  elements  of  the  PEN  and 
the  AST  can  be  obtained  by  further  differentiations: 
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where  Xi=r2,  yi=0,  dz2  =  0.  For  the  arc  ^12, 
x2  =  rcos0 ,  y2  =  rsm0,  dy2  =  r cos OdO,  0  e[0m,  -0m],  in  which 
0m  =  cos-1(ri/r2)  +  cos-1(ri/r).  For  the  line  23,  x2  =  rcosOm, 
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Fig.  4.  Configuration  of  two  sets  of  equally  spaced  curve  segments  extending  along  two  parallel  generating  lines. 
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For  the  non-concave  exterior  surface  of  the  AST  and  the  PEN, 
dPdASTi,dAST2,out  =  0»  an4  dFdPENi.dPEN2.out =  0-  With  the  reciprocity 
rule,  dFdAST1  0Ut  dpEN2/dFdPEN1  dAST2i0ut =  ^- 

The  view  factor  between  the  ring  elements  on  the  interior  sur¬ 
face  of  the  AST  (dFdAST1  dAST2  in)  and  the  view  factor  of  the  ring 
element  to  the  end  disk  (FdAST  in  end)  can  be  obtained  by  differenti¬ 
ating  the  view  factor  of  the  disk  to  the  parallel  coaxial  disk. 

Define  Z=z/2/rASTin,X=  |Z2-Z!|,  then 


dFdASTl,  dAST2,  in 


_  X(2X2  +  3) 
2(X2  +  1  )3/2 
2X2  +  1  x 
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dX2 ,  fdAST  ,  in,  end 


(15) 


3.  Analytical  finite-difference  view  factors 

The  differential  view  factor  model  generally  requires  differen¬ 
tial  elements  small  enough  to  obtain  an  accurate  calculation.  For  a 
typical  SOFC,  the  cell  length  is  usually  much  bigger  than  its  char¬ 
acteristic  size,  i.e.  L»Dch  and  Wch  for  PSOFC;  and  L»rou t  in 
TSOFC.  So  a  considerable  number  of  grids  are  required  in  distributed 
SOFC  modeling,  which  leads  to  expensive  computation.  On  the 
other  hand,  there  are  numerical  discontinuities  in  the  differential 
view  factors  of  PSOFC  (e.g.  /i(X,Y)^  oo  as  Y=0  in  Eq.  (4)),  which 
collapses  the  overall  heat  balance  and  generates  unreasonable  tem¬ 
perature  distributions.  The  analytical  finite-difference  view  factor 
is  a  good  alternative  to  these  problems. 

Fig.  4  shows  a  configuration  of  two  sets  of  equally  spaced  curve 
segments  extending  along  two  parallel  generating  lines.  According 
to  the  reciprocity  and  additivity  rules,  the  view  factor  from  the  con¬ 
tinuous  segments  a2  ~  ai+k  to  the  continuous  segments  b[  ~  bi+k  can 
be  obtained  by 

1  (  I+k-1 

Farai+k,brbi+k  =  I  FFarai+k_A,brbi+k_A  +  ^  Fai+k,bj 

V  j=i 


Because  of  symmetry,  the  view  factor  from  the  element  a2  to  the 
element  bj  is  only  dominated  by  their  relative  positions.  Define 

Fahbj  =  Fai+k,bj+k  =s(\i~j\)  >  Fanai+k,bnbi+k  =fU<  +  1)  (17) 

and  substitute  it  into  Eq.  (16)  to  generate: 

k 

(k  +  l)/(k  +  l)  =  k/(k)  +  2^g(i)+/(l)  (18) 

1=1 
So, 

g(o)=nn  g(i)=/(2)-/(i) 

S(k)  =  l  [(k  + 1  ]f(k  + 1 )  -  2 k/(k)  +  (k  -  1  )/(k  -  1 )]  (19) 

(k  =  2 ...  n  -  1) 


3.1.  PSOFC  configuration 


Considering  a  rectangular  cuboids  enclosed  by  the  first  k  PEN 
elements  PEN(l)~PEN(/<)  (denoted  by  Ai),  the  directly  opposed 
CON  elements  CON(l)~CON(k)  (denoted  by  A2)  and  the  two  bot¬ 
toms  (denoted  byA3  andA4;A3  corresponding  to  the  first  element), 
the  view  factors  F13(/<)  and  F34(/<)  can  be  calculated  using  the 
equations  for  C-l  1  and  C-14  configurations  in  the  view  factor  cata¬ 
logue  [11].  The  following  equation  can  be  generated  when  defining 
X=Dch/Wch,  Al  =  L/n,  Y=AllWch: 


/pen.conM  =  1  -  2F1?3(/<)  (/<  =  1 ...  n) 


/con,  con  (k)  = 
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(20) 


The  view  factors  FPEN(i),coN(j)  an4  ^con(i).cono')  PSOFC  can  then 
be  calculated  by  substituting  Eq.  (20)  into  Eq.  (19).  The  view  factors 
from  the  kth  PEN  element  or  CON  element  to  the  top  end  (top)  and 
bottom  end  (bot)  of  the  groove  can  be  obtained  by 

PpEN,top(l  )  —  Pi, 3(1 ) 

Fpen,  top  (/<)  =  kFh3(k)  -  (k  -  1  )Fh3(k  -  1)  (k  =  2...n) 
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PpEN,  bot(^)  =  PPEN,  top(n  +  1  -  k) 

PcoN,bot(/<)  =  PcoN,top(n  +  l  -k)  {k  =  1  ...n) 

3.2.  TSOFC  configuration 
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icTT 


i+k 


j=i 


Considering  a  configuration  of  coaxial  cylinders  enclosed 
(16)  by  the  first  k  AST  exterior  elements  AST,out(l)~AST,out(/<) 
(denoted  by  Ai),  the  directly  opposed  PEN  interior  elements 
PEN,  in(l)~PEN,  in(/<)  (denoted  by  A2)  and  the  two  annular  disks 
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Fig.  5.  Axial  distribution  of  the  dimensionless  tube  temperature  among  numerical  solutions  based  on  the  differential  view  factor  (VF)  model,  Siegel’s  and  Perlmutter’s  method 
and  the  finite-difference  view  factor  model  at  different  ratio  of  tube  length  to  tube  diameter  L/D  and  different  grid  number  n. 


at  base  or  top  (denoted  by  A3  and  A4 ,  and  A3  is  correspond¬ 
ing  to  the  first  element),  the  view  factors  F2> 2(/<),  F2ii(/<),  F23(/<) 
and  F1>3(/<)  can  be  respectively  calculated  using  the  equations  for 
C-91,  C-92,  C-93  and  C-77  configurations  in  the  view  factor  cat- 
alogue  [11],  The  view  factors  FpEN,in(f),PEN>in[/)  and  FpEN,in(i),AST,out(j) 
in  TSOFC  can  then  be  obtained  by  substituting /PENipENiin(F)  =  F2>2(/<) 
and/PENiAST,out(^)  =  F2,i(/<)  into  Eq.  (19).  The  view  factors  from  the 
kth  AST  exterior  element  or  PEN  interior  element  to  the  top  and 
bottom  annular  ends  can  be  calculated  by 

FpEN,  top(l  )  —  F2?3(l  ) 

FpEN,top(/<)  =  /<F2?3(/<)  —  (/<  —  1)F2j3(/<  -  1)  ( k  =  2...n )  v  j 

Fast,  out,  top(^ )  — Fi?3(1  ) 

Fast,  out,  toP(/<)  =  fcFif3(*)  -  (k  -  1  )Fh3(k  -  1)  (k  =  2  ...  n)  lZDj 

FpEN,  bot(F)  =  Fpen,  top(n  +  1  -l<) 

Fast,  out,  bot(F)  —  Fast,  out,  top(n  +  !—/<)  (k  =  1 ...  n) 


for  the  gas  and  solid  phases  are,  respectively  [12,13]: 

PgUgcp,g-^  —  -  Fg),  Fg|z=o  =  Tg  m  (28) 

0  =  -^  (vT?  -Bs)+h(Tg-Ts)  +  q  (29) 

where  z  is  the  axial  length  coordinate  measured  from  the  tube 
entrance,  Fg  is  the  gas  temperature,  Fg  in  is  the  inlet  gas  temper¬ 
ature,  pg  is  the  gas  density,  ug  is  the  mean  gas  velocity,  cp>g  is 
the  mass  specific  heat  of  gas,  h  is  the  convective  heat-transfer 
coefficient,  D  is  the  tube  diameter,  Ts  is  the  tube  temperature, 
<7  =  5.67  x  1 0-8  W  nrr2 1<-4  is  the  Stefan-Boltzmann  constant,  s  is 
the  emissivity  of  solid  phase,  q  is  the  heat  added  per  unit  area  at 
the  tube  wall.  Note  that  the  conductive  heat  transfer  is  not  included 
in  Eq.  (29).  The  total  outgoing  radiation  per  unit  area  from  a  tube 
surface  element  Bs  is  related  to  the  tube  temperature  and  view 
factors: 


The  view  factor /AST,AST,in(^)  from  the  continuous  AST  interior 
elements  AST,  in(  1 )  ~  AST,  in(/<)  to  itself  can  be  calculated  using  the 
equation  for  the  C-78  configuration  in  the  view  factor  catalogue 
[11].  The  view  factor  FAsT,in(i),AST,in(/)  can  then  be  calculated  by  Eq. 
(19)  and  the  view  factor  from  the  kth  AST  interior  element  to  the 
circular  end  can  be  obtained  by 

Fast,  in,  end  ( 1 )  =  0  ^ •  5  1  -  /asT,  AST,  in  ( 1 )] 

Fast,  in,  end (F)  =  0.5  1  -  /</ast,  AST,  in(F)  +  (F  -  1  Vast,  AST,  in(F  -  1 )] 

(/<  =  2 ...  n) 

(27) 


4.  Simulation  and  discussions 

4.1.  Problem  of  convective  and  radiant  heat  transfer 

For  problems  of  convective  and  radiant  heat  transfer  for  flow  of  a 
transparent  gas  in  a  tube  with  a  gray  wall,  the  governing  equations 


Bs(z)  =  £oTs4  +  (l-e) 


[  Bs(?)dFz,j  +  are4inF2jin  +  are4  outFz, 

Jo 


out 

(30) 


where  dFz  ^  =  I<(  \z  -  §|)d§  is  the  differential  view  factor  between  two 
ring  elements  at  location  z  and  §,  I<(\z-$\)  is  the  geometric  kernel, 
Fz  in  and  FZi0Ut  are  the  view  factors  from  an  element  at  location  z  to 
the  circular  opening  at  the  inlet  and  outlet  ends  of  the  tube,  respec¬ 
tively.  For  a  circular  tube,  I<(\z  -  §|),  Fz  in  and  FZi0Ut  can  be  calculated 
using  Eq.  (15).  Tein  and  re,out  are  the  effective  black-body  tem¬ 
perature  of  the  environment  at  the  inlet  and  outlet  end  openings, 
which  are  usually  set  to  be  the  inlet  and  outlet  gas  temperatures, 
he.  Fe  in  =  Fg  in,  Fe ^  =  Fg)0ut =  Fg  |z=l  [  13  ]. 

An  additional  check  for  each  numerical  solution  can  be  made 
by  making  sure  that  it  always  satisfies  an  overall  heat  balance.  The 
total  input  heat  to  the  tube  wall  Qjn  is 


Qin  —  7tD(qL  + 


/  <infz,in dz+  [ 

Jo  Jo 


crT; 


e,  outrz,  ou 


t  dz) 


(31) 
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Fig.  6.  Comparison  between  the  finite-difference  (FD)  and  first-order  approximation  of  differential  (Diff)  view  factors  (a)  between  interior  tube  ring  elements  (b)  from  an 
interior  tube  ring  element  to  the  circular  end  at  different  grid  number  n.  The  relative  element  number  equal  to  zero  means  two  directly  opposed  elements  and  the  element 
sequence  number  equal  to  1  represents  the  first  element  adjacent  to  the  inlet  opening  end. 


and  the  total  output  heat  from  the  tube  wall  Qout  is 
Qout  =  ^7rD2pgUgCp?  g(Tg;  out  —  Tg  in) 


+  7 ID 


j  Bs(z~)Fz  jndz  +  J  Bs(z)Fz,  outdz^ 


(32) 


The  overall  heat  balance  is  satisfied  when  Qin  =  Qout;  the  relative 
error  on  the  overall  heat  balance  is  then  defined  as  |1  —  Qout/Qinl- 

When  defining  the  dimensionless  coordinate  as  Z=z/D,  Stan¬ 
ton  number  as  S  =  4/i/(pgUgCP)g),  dimensionless  temperature  as 
t=T(or/q)1/4  and  the  dimensionless  heat-transfer  coefficient  as 
H  =  (hlq)(qlcry14,  Eqs.  (28)-(32)  can  be  reduced  to  dimensionless 
form  [13].  Siegel  and  Perlmutter  further  transformed  the  integral 
equation  (30)  to  a  second-order  differential  equation  by  use  of  an 
approximate  exponential  function  for  the  geometric  kernel  I<  [13]. 
However,  it  requires  complex  mathematical  derivations  to  find  the 
two  boundary  conditions.  When  the  conductive  heat  transfer  is 
considered  in  Eq.  (29),  it  will  generate  a  fourth-order  equation. 

Fig.  5  shows  the  axial  distribution  of  the  dimensionless  tube 
temperature  for  three  different  numerical  solutions.  The  first  solu¬ 
tion  is  from  the  dimensionless  form  of  Eqs.  (28)-(30)  with  the 
differential  view  factors  in  Eq.  (15),  the  second  solution  is  from 
the  second-order  differential  system  of  Siegel  and  Perlmutter, 
and  the  third  solution  is  obtained  from  manual  discretization  of 
Eqs.  (28)-(30)  with  the  finite-difference  view  factors  in  Eqs.  (19) 
and  (27).  With  the  following  parameters  s  =  0.8,  q  =  1  x  1 04  W m-2, 
S  =  0.00988,  H  =  0.8,  tgin  =  1.5,  all  three  numerical  solutions  were 
calculated  using  the  commercial  software  gPROMS,  which  pro¬ 
vides  equation-based  solvers  for  differential  and  integral  hybrid 
systems  [14].  As  shown  in  Fig.  5(a),  for  a  short  tube  with  L/D  =  5, 
all  three  numerical  methods  led  to  an  identical  temperature  distri¬ 
butions  at  a  default  grid  number  n  =  50,  and  the  relative  errors  on 
the  overall  heat  balance  were  0.0245,  5.24  x  10-4  and  5.55  x  10-16 
for  the  differential  view  factors  solution,  Siegel  and  Perlmutter 
solution,  and  finite-difference  view  factors  solution,  respectively. 
When  L/D  =  50  and  n  =  50,  the  element  size  is  comparable  to  the 
tube  diameter  dz  =  L/n  =  D,  but  it  then  becomes  difficult  to  obtain 
a  solution  that  converges  when  using  the  differential  view  factors 


method.  As  shown  in  Fig.  5(b),  the  temperature  distribution  based 
on  finite-difference  view  factors  is  smoother  than  that  obtained 
using  Siegel’s  and  Perlmutter’s  computation,  while  the  relative 
error  on  the  overall  heat  balance  2.49  x  1 0-13  is  much  smaller  than 
that  for  the  Siegel’s  and  Perlmutter  method  (0.00169).  In  the  com¬ 
mercial  Westinghouse  TSOFC  design,  the  size  of  the  air  supply  tube 
is  Fast, in  =  2.5  x  10-3  m  and  its  length  is  L  =  1.5m,  which  means  a 
high  L/D  ratio  of  300.  In  this  case,  a  numerical  convergence  using 
Siegel’s  and  Perlmutter’s  method  also  requires  a  large  grid  number 
to  get  a  small  element  size.  When  n  =  150  and  n  =  200,  the  relative 
errors  on  the  overall  heat  balance  using  the  Siegel’s  and  Perlmut¬ 
ter’s  method  are  up  to  0.538  and  0.164,  respectively.  As  shown  in 
Fig.  5(c),  the  computation  based  on  the  finite-difference  view  fac¬ 
tors  displays  a  different  temperature  distribution  for  the  second 
half  of  tube.  It  took  only  30  ms  for  the  computation  based  on  the 
finite-difference  view  factor  model  to  find  a  solution  with  a  very 
small  relative  error  on  the  overall  heat  balance  1.177  x  10-14  by 
keeping  a  default  grid  number  n  =  50  at  L/D  =  300. 

The  above  problem  can  be  further  explained  when  comparing 
the  finite-difference  and  differential  view  factors  methods.  Accord¬ 
ing  to  the  physical  meaning  of  the  view  factors,  their  value  should 
be  not  greater  than  unity.  As  shown  in  Fig.  6(a),  the  finite-difference 
view  factors  between  the  interior  AST  elements  always  follow  this 
rule,  while  the  differential  view  factors  appear  to  be  always  greater 
than  1  for  the  two  directly  opposed  elements  even  for  n  =  200.  Here, 
the  differential  view  factors  are  calculated  by  approximation  of  the 
first-order  differences,  i.e.  dFz  ^  =  I<{  |z  -  §|  )d§  ^  KL/n.  For  large  value 
of  L/D  and  small  grid  number  n,  the  element  size  is  comparable 
to  the  characteristic  size  (dz  =  L/n^D),  which  leads  to  computa¬ 
tional  problem  with  differential  view  factors.  As  shown  in  Fig.  6(b), 
the  view  factor  from  the  first  infinitesimal  interior  AST  element  to 
the  circular  end  is  equal  to  0.5,  independently  of  the  grid  num¬ 
bers  (FdAST.in.end  =  0.5  as  X=0  in  Eq.  (15)),  and  the  corresponding 
finite-difference  view  factors  are  much  lower  than  this  limiting 
value  (FAST,in,end(l)  =  0-5  as  oo).  Fig.  7  shows  similar  results  for 
the  differential  and  finite-difference  view  factors  for  TSOFC  with 
t ast, out  =  4x1 0-3  m,  rin  =  8.66  x  1 0-3  m  and  L  =  1 .5  m. 

Numerical  problems  are  more  pronounced  in  the  computation 
of  radiation  heat  transfer  in  planar  SOFCs,  because  of  discontinu- 
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Fig.  7.  Comparison  between  the  finite-difference  (FD)  and  first-order  approximation  of  differential  (Diff)  view  factors  (VF)  (a)  from  PEN  element  to  PEN  element,  (b)  from 
PEN  element  to  AST  element,  (c)  from  PEN  element  to  the  annular  end  and  (d)  from  AST  element  to  the  annular  end  in  the  tubular  SOFC  configuration.  The  relative  element 
number  equal  to  zero  means  two  directly  opposed  elements  and  the  element  sequence  number  equal  to  1  represents  the  first  element  adjacent  to  the  inlet  opening  end. 


ities  in  the  geometric  kernel  when  considering  differential  view 
factors.  With  a  typical  PSOFC  geometry  [15],  Wch  =  3xl0-3m, 
Dch  =  1  x  10-3  m,  L  =  0.1  m,  Fig.  8  shows  a  comparison  between  the 
finite-difference  view  factors  and  the  first-order  approximation  of 


the  differential  view  factors  in  the  PSOFC  configuration.  As  shown  in 
Fig.  8(a)  and  (b),  the  finite-difference  view  factors  between  PEN  and 
CON  elements  are  always  lower  than  unity,  while  the  differential 
view  factors  between  two  directly  opposed  elements  are  infinite 
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Fig.  9.  Distribution  of  the  current  density  and  solid  temperature  from  the  distributed  cell  model  with  interface  (INT-)  type  PEN  model  (assuming  that  electrochemical  reactions 
only  occur  at  the  electrode/electrolyte  interface,  and  reforming  and  water  gas  shift  chemical  reactions  only  occur  at  the  anode/flow  channel  interface)  for  90%  H2  -1 0%  H20  fuel 
system,  (a)  and  (b)  under  co-flow  condition,  (c)  and  (d)  under  counter-flow  condition.  The  geometries  of  PSOFC  are:  cathode  thickness  <5a  =  8C  =  50  |jim,  electrolyte  thickness 
8e  =  1 50  |jim,  channel  width  Wch  =  3  mm,  channel  depth  Dch  a  =  Dch  c  =  1  mm,  rib  width  Wrib  =  2.42  mm,  total  bipolar  plate  depth  Dt,a  =  Dt>c  =  2.5  mm,  channel  length  1  =  100  mm, 
and  the  number  of  channels  nch  =  18  [15]. 


Cfi(X,Y)^oo  and/2(X,y)^  oo  as  Y=0  in  Eqs.  (4)  and  (5)).  As  shown 
in  Fig.  8(c)  and  (d),  the  view  factor  from  the  first  infinitesimal  PEN 
and  CON  elements  to  the  opening  end  is  equal  to  0.5  at  different  grid 
numbers  (FdPEN,end  =  0-5  and  FdcoN.end  =  0-5  as  y=0  in  Eqs.  (6)  and 
(7)),  and  the  corresponding  finite-difference  view  factors  are  much 
lower  than  this  limiting  value  (FpEN,end(l ) =  0-5  and  FcoN,end(l ) =  0-5 
as  n  oc). 


4.2.  Application  to  distributed  SOFC  model 


The  energy  conversation  in  the  PEN  is  [10]: 

9FPen  f/.  d2rPEN  c  at  i  u  (T  ■ 

PPENCp,  pen  df  =  kPEN— 7^2 - ^m,  PEN,  c  /  Nj ,c|C/Crfc j(lc. 


(n,-,; 


-  N, 


uL/cIHiTa) 

tfa,i(’FpEN)J  -  Sm,  pen,  cK^)V cell 


+  5h’ PEN>  k^k’  pen(^i<  -  Fpen)  -  5h>  PEN’  /<^pen,  k  (33) 

k= a,  c  k= a,  c 


where  Ppen.  cp, pen,  ^pen  are  the  density,  mass  specific  heat  capacity, 
and  heat  conductivity  of  the  PEN;  Sm  PEN,/<  and  5h  PEN,k  are  the  effec¬ 
tive  mass  and  heat  transfer  area  per  unit  volume  of  PEN  between 
gas  and  solid  phases  in  the  anode  (/<  =  a)  and  cathode  (/<  =  c)  flow 
channels;  Ta,  Tc,  TPEn  are  the  temperatures  of  fuel,  air  and  PEN,  h 
is  the  corresponding  convective  heat  transfer  coefficient,  Hj  is  the 
molar  specific  enthalpy  of  species  i,  Nz  is  the  molar  flux  of  species 
i  at  the  electrode/flow  channel  interface  (anode  side:  A/C;  cathode 
side:  C/C),  /( z)  is  the  local  current  density,  and  Vcell  is  the  cell  voltage 
which  is  assumed  uniformly  distributed  along  the  z  coordinate. 

On  the  right  hand  side  of  Eq.  (33),  the  first  item  is  heat  con¬ 
duction,  summations  of  the  second  and  third  items  are  the  total 
electrochemical  and  chemical  reaction  heats,  the  fourth  item  is 
the  electric  power,  the  fifth  item  is  the  convective  heat  transfer 


between  the  gas  phase  and  the  PEN,  and  the  sixth  item  is  the  radi¬ 
ation  heat  transfer.  For  gray  bodies,  the  radiation  heat  flux  out  of 
the  PEN  element  surface,  qPEN  is  obtained  by 


4pen ,k(z) 


£pen  [orF^EN(z)-Bp£N(2:)] 

- - — - -  (/<  =  a,  c) 

1  -  £PEN 


(34) 


For  PSOFC  configuration, 


FpEN ,kiz)  =  ^PENOrFp£N(z)  +  (1  -  G  PEN)tfFeFdPENz,end,/< 

/  FpEN(*)dFdpENz,  dPENx,  k 

Jo 


+  (1  -  £pen) 
rL 


JO 


FcON,  k(*)dFdPENz,  dCONx,  k 


(35) 


FcON,  k(z)  =  £C0N^F^0N(z)  +  (1  -  £con)  [^FgFdcoNz,  end,  k 
+  /  FCOn,  k(x)dF dCONz,  dCONx,  k 

Jo 

+  /  Fpen,  k(*)dFdC0Nz,  dPENx,  k 

Jo 


(36) 


where  £pen  and  £Con  are  the  emissivities  of  PEN  and  CON,  BPEn 
and  Bcon  are  the  total  outgoing  radiation  per  unit  area  from  the 
PEN  and  CON  element  surfaces,  respectively,  and  Te  is  the  effec¬ 
tive  black-body  temperature  of  the  environment,  which  is  set  to  be 
the  gas  temperature  at  the  corresponding  end  of  the  flow  channel 
[13].  Note  that  the  analytical  differential  view  factors  are  used  here, 
which  should  be  replaced  by  the  analytical  finite-difference  view 
factors  when  using  the  discrete  governing  equations.  The  energy 
balances  for  the  CON  in  PSOFC  and  AST  in  TSOFC  are  simpler  and 
can  be  obtained  similarly  [10]. 

In  contrast  to  the  detailed  radiation  model  mentioned  above, 
only  surface-to-surface  radiation  exchange  between  two  directly 
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Fig.  10.  Distribution  of  the  current  density  and  solid  temperature  from  the  distributed  cell  model  with  INT-type  PEN  model  and  control  volume  (CV-)  type  PEN  model 
(electrochemical  and  chemical  reactions  occur  everywhere  in  electrodes)  for  30%  pre-reformed  methane-fueled  system,  (a)  and  (b)  under  co-flow  condition,  (c)  and  (d)  under 
counter-flow  condition.  The  geometries  of  PSOFC  are  the  same  as  those  in  Fig.  9. 


opposed  elements  is  taken  into  account  in  most  of  the  state-of- 
the-art  SOFC  models  [4-6].  For  an  assumed  enclosed  system  of  PEN- 
CON  rectangle,  the  simple  radiation  model  means 

a  = _ (rPEN  ~  Tcon) _  (k  =  a  cl 

PEN'  k  l/ePEN  +  Wch(l/£C0N  -  l)/(Wch  +  2Dch  j() 

(37) 

Our  SOFC  model  was  validated  by  comparison  with  experimen¬ 
tal  data  from  the  IEA  Benchmark  Test  for  a  typical  PSOFC  design 
[15].  For  both  FI2-FI2O  and  pre-reformed  methane  fuel  systems, 
the  operating  pressure  is  1  bar,  the  air  and  fuel  inlet  temperatures 
are  900  °C,  the  air  stoichiometry  is  7,  and  the  overall  fuel  utilization 
is  85%.  When  the  analytical  differential  view  factors  are  used,  the 
built-in  discrete  methods  of  gPROMS  fail  to  obtain  a  solution  with 
a  satisfactory  relative  error  on  the  overall  heat  balance  because  of 
the  discontinuity  problem  in  the  PSOFC  configuration.  Based  on  the 
analytical  model  of  finite-difference  view  factors  and  stagger  grids, 
a  fast  (average  10  s  for  50  discrete  grids)  and  accurate  solution  is 
obtained  with  a  relative  error  of  the  overall  heat  balance  <0.1  %  [  1 0]. 

Fig.  9  shows  the  distribution  of  the  current  density  and  PEN 
temperatures  under  both  co-flow  and  counter-flow  conditions  for 
90%H2-1 0%Fl2O  fuel  system.  Under  the  co-flow  condition,  as  shown 
in  Fig.  9(b),  the  solid  temperature  predicted  based  on  the  “detailed” 
radiation  model  is  almost  10  K  lower  than  the  “simple”  calculation. 
More  importantly,  the  smaller  solid  temperature  gradient  implies  a 
smaller  thermal  stress  from  the  “detailed”  prediction  than  from  the 
“simple”  prediction.  Under  the  counter-flow  condition,  as  shown  in 
Fig.  9(d),  there  is  an  obvious  nonlinear  distribution  of  the  PEN  tem¬ 
perature  in  the  “detailed”  calculation,  which  is  not  well  captured  in 
the  “simple”  calculation.  Because  of  the  strong  radiant  heat  loss  to 
the  environment  in  the  “detailed”  calculation,  the  maximum  neg¬ 
ative  gradient  of  the  solid  temperature  appears  close  to  the  fuel 
inlet.  Compared  to  the  “simple”  calculation,  the  “detailed”  calcu¬ 
lation  avoids  the  overestimation  of  the  maximum  current  density 
and  solid  temperature  (especially  under  counter-flow  condition), 


which  provides  a  more  uniform  distribution  of  the  current  density 
as  shown  in  Fig.  9(a)  and  (c). 

Fig.  10  shows  the  distribution  of  the  current  density  and 
solid  temperature  under  co-flow  and  counter-flow  conditions 
for  fuel  system  of  1 7.1  %CH4-26.26H2-49.34%H20-2.94%CO-4.36% 
C02,  which  means  30%  pre-reformed  gas  at  CO  shift  equilibrium 
resulting  from  a  methane/water  mixture  with  a  molar  ratio  of 
FI2O/CFI4  =  2.5.  Similar  to  the  results  with  the  FI2-FI2O  fuel  system, 
the  distributions  of  the  current  density  and  solid  temperature  from 
the  “detailed”  calculation  are  more  uniform  than  those  from  the 
“simple”  calculation.  The  maximum  PEN  temperature  and  cell  per¬ 
formance  in  the  “detailed”  calculation  are  also  generally  lower  than 
in  the  “simple”  calculation.  Unlike  the  results  in  Fig.  10(b)  and  (d) 
shows  a  non-monotonic  distribution  of  PEN  temperature  in  both 
the  “simple”  and  “detailed”  calculation  under  counter-flow  con¬ 
ditions  because  of  the  strongly  endothermic  methane  reforming 
reaction.  By  considering  the  radiant  heat  loss  to  the  environment, 
the  “detailed”  radiant  model  led  to  a  more  intense  negative  gra¬ 
dient  in  solid  temperature  close  to  the  fuel  inlet.  As  a  result, 
Fig.  10(c)  shows  a  nonlinear  distribution  of  the  current  density  in 
the  “detailed”  calculation. 

5.  Conclusions 

Radiation  heat  transfer  plays  an  important  role  in  temperature 
distribution  in  solid  oxide  fuel  cells  (SOFC).  A  SOFC  model  with  a 
detailed  radiation  model  is  necessary  for  accurately  predicting  the 
cell  performance.  For  the  two  typical  designs  of  tubular  and  planar 
configurations,  a  mathematical  model  of  view  factors  for  radiation 
heat  exchange  in  longitudinally  distributed  SOFC  modeling  was 
introduced  in  this  paper. 

An  analytical  model  of  differential  view  factors  was  first  devel¬ 
oped  for  the  configurations  of  diffusive  interchange  in  planar  and 
tubular  SOFCs.  An  accurate  solution  is  usually  obtainable  when  the 
size  of  the  differential  elements  is  small  enough.  Flowever,  the 
cell  length  is  generally  much  larger  than  its  characteristic  size  in 
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a  typical  SOFC.  So  a  large  amount  of  discrete  grids  are  required 
to  ensure  numerical  accuracy,  which  leads  to  expensive  computa¬ 
tion.  In  the  case  of  planar  configuration,  the  differential  view  factor 
model  failed  because  of  the  discontinuities  of  some  geometric  ker¬ 
nels.  An  analytical  model  of  finite-difference  view  factors  was  then 
developed  to  resolve  the  non-convergence  problem  of  the  solution. 

The  mathematical  model  of  view  factors  was  first  validated  via 
comparison  of  known  results  of  a  classical  problem  of  convective 
and  radiant  heat  transfer  in  the  case  of  flow  of  a  transparent  gas  in 
a  tube  with  a  gray  wall.  All  computational  results  were  compared 
to  Siegel’s  and  Perlmutter’s  solution  by  transforming  the  integral 
equation  of  radiant  energy  balance  into  a  second-order  differential 
equation  [13].  For  a  small  ratio  of  tube  length  to  tube  diameter  (L/D), 
the  differential  view  factor  model  leads  to  an  accurate  solution.  For 
a  large  L/D  ratio,  numerical  convergence  using  Siegel’s  and  Perlmut¬ 
ter’s  method  requires  a  large  grid  number  to  get  a  small  element 
size.  A  fast  and  accurate  computation  is  available  for  the  finite- 
difference  view  factor  model  with  a  very  small  relative  error  for 
the  overall  heat  balance  by  keeping  a  default  grid  number  (n  =  50) 
without  a  limitation  of  large  L/D  ratio.  Another  advantage  of  the 
finite-difference  view  factor  model  is  that  it  avoids  extra  complex 
mathematical  derivations  of  the  governing  equations  like  the  ones 
Siegel  and  Perlmutter  did. 

In  addition  to  a  simple  radiation  model,  a  detailed  radiation 
model  based  on  the  analytical  view  factors  was  then  applied  to  an 
in-house  distributed  cell  model  and  solved  using  the  commercial 
software  gPROMS  for  a  typical  PSOFC  configuration  [10].  Com¬ 
pared  to  the  simple  radiation  model,  which  only  takes  into  account 
the  surface-to-surface  radiation  exchange  between  two  directly 
opposed  elements,  the  distributed  cell  model  with  the  detailed 
radiation  model  predicts  more  uniform  distributions  of  the  local 


current  density  and  solid  temperature,  which  also  means  a  smaller 
thermal  stress.  In  principle,  the  analytical  view  factors  models  can 
be  applied  to  other  types  of  cell  designs  with  regular  geometry, 
such  as  the  flattened  tubular  SOFC  and  integrated  planar  SOFC. 
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